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Anisotropic mesh adaptation for 3D anisotropic 
diffusion problems with application to fractured 

reservoir simulation* 

Xianping Li^ Weizhang Huang^ 


Anisotropic mesh adaptation is studied for linear finite element solution of 3D anisotropic 
diffusion problems. The M-uniform mesh approach is used, where an anisotropic adaptive 
mesh is generated as a uniform one in the metric specified by a tensor. In addition to mesh 
adaptation, preservation of the maximum principle is also studied. Some new sufficient 
conditions for maximum principle preservation are developed, and a mesh quality measure 
is defined to server as a good indicator. Four different metric tensors are investigated: 
one is the identity matrix, one focuses on minimizing an error bound, another one on 
preservation of the maximum principle, while the fourth combines both. Numerical ex¬ 
amples show that these metric tensors serve their purposes. Particularly, the fourth leads 
to meshes that improve the satisfaction of the maximum principle by the finite element 
solution while concentrating elements in regions where the error is large. Application of 
the anisotropic mesh adaptation to fractured reservoir simulation in petroleum engineering 
is also investigated, where unphysical solutions can occur and mesh adaptation can help 
improving the satisfaction of the maximum principle. 
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1 Introduction 

We are concerned with the linear finite element solution of the three dimensional boundary value 
problem (BVP) of the diffusion equation 

-V-(DVix) = /, in n (1) 
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subject to the Dirichlet boundary condition 

u — on dO. (2) 

where 12 is a bounded polyhedral domain, / and g are given functions, and D is the diffusion matrix. 
We assume that D = D(x) is a symmetric and uniformly positive definite matrix-valued function on 
12. It includes both isotropic and anisotropic diffusion as special examples. In the former case, D takes 
the form D = a{x)I^ where / is the 3x3 identity matrix and a — a{x) is a scalar function. In the 
latter case, on the other hand, D has not-all-equal eigenvalues at least on a portion of 12. 

Anisotropic diffusion problems arise from various branches of science and engineering, including 
plasma physics [211 ESI E3 ED ET], petroleum engineering ID El ED ED], and image processing uni 
[ID EH EH ED- When a conventional numerical method is used to solve the problems, spurious 
oscillations may occur in computed solutions. Numerous research has been done for two dimensional 
(2D) problems; among other works, we mention a few here, [HI [121 [El ESI ESI ESI ESI SOI El El El Bl 
ESI EH [Ml ESI [66]. A common approach is to design a proper discretization method and/or a proper 
mesh so that the numerical solution satisfies the maximum principle (MP). Recently, an anisotropic 
non-obtuse angle condition was developed in ISDIinillDllSl for the linear finite element solution of 
both time independent and time dependent anisotropic diffusion problems to satisfy MP. 

On the other hand, much less work has been done for 3D anisotropic diffusion problems. Although 
MP preservation has been studied in general dimensions e.g. in |El[I2lE3IMlEHlEaiinillD , most 
of them either consider isotropic diffusion or present numerical examples only in ID and 2D. For 
example, only isotropic diffusion is considered in II21EZ!. It is shown in [39| that the 3D Delaunay 
triangulation does not generally produce a mesh with which the numerical solution satisfies MP. Mesh 
conditions are studied in [8] for a reaction-isotropic-diffusion problem for general dimensions and 
numerical examples in ID and 2D are presented. The difficulty of MP satisfaction for 3D problems is 
remarked in both [8| and [39] . 

The objective of this paper is to study the linear finite element solution of 3D anisotropic diffusion 
problems. The focus will be on MP preservation and mesh adaptation. Four different metric tensors 
used in anisotropic mesh generation will be considered. This study is a 3D extension of the work 
m- Moreover, new sufficient conditions will be developed for the linear finite element approximation 
to satisfies MP, and a mesh quality measure will be defined to provide a useful indication for MP 
satisfaction. The mesh quality measure is developed along the approach of [28]. But the interested 
reader is also referred to [351 |3B] and references therein for different mesh optimization methods 
and quality metrics. Furthermore, the application to fractured reservoir simulation in petroleum 
engineering will also be investigated, where unphysical solutions can occur and mesh adaptation can 
help improving MP satisfaction. 

An outline of the paper is given as follows. The linear finite element formulation for BVP 0 and 
Q is given in Section]^ MP preservation and some sufficient conditions will be discussed. Section]^ 
contains the discussion of anisotropic mesh adaptation and four metric tensors. Numerical examples 
are given in Section followed by the investigation of application of anisotropic mesh adaptation to 
fractured reservoir simulation in Section El Conclusions are drawn in Section El 

2 Finite Element Formulation 

In this section we brieffy describe the piecewise linear finite element discretization for BVP Q and 
0 and state a few properties of the discretization. 
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Let Ug — {v ^ I v\qq = g}. The weak formulation of BVP Q and Q is to find u ^Ug such 

that 

/ (Vv)'^I])'Vudx = / fvdx^ \/v^Uo. 

Jq Jq 

For the finite element discretization, we assume that a tetrahedral mesh Th is given on Q. Let 
be the piecewise linear interpolation of g on the boundary dQ and be the piecewise linear finite 

element space on Th with the boundary data g^. Then, the finite element formulation for BVP Q 
and (2) is to find G 


such that 


/ (V^ 

Jq 


MT 


B\/u^dx= [ fv^dx, 

Jq 


G U^. 


(3) 


This discretization is standard and it is expected that converges to 1 / at a rate of second order in 
the norm and first order in norm. We take the reference element K to be a unitary equilateral 
tetrahedron and denote an element in the mesh Th by K. Moreover, denote by {i = 1, ...,4) the 
unit inward normal to the face facing the vertex of K. Let Fk be the affine mapping from K to 
K and F'k the Jacobian matrix of Fk- 

Theorem 2.1 (Li and Huang uni). If the mesh satisfies 


< 0, Vi ^ j, i,j = ViV € % (4) 

where is the average of D over K, then the linear finite element seheme for solving BVP ^ 
and satisfies the maximum prineiple, 


f < 0 in 


max u^ix) 
xeQudQ 


max u^ix). 
xedQ 


Preserving MP is crucial to avoiding artificial oscillations in the computed solution. It is thus 
interesting to know what meshes satisfy the anisotropie non-obtuse angle eondition Obviously, 
a sufficient condition is 


{F'K)-^JiK{F'K)-^ = CkI, VK € Th 


(5) 


where ck is a positive scalar constant on K and / is the 3x3 identity matrix. 

A weaker condition is stated in the following theorem, where a general case for the d space dimension 
{d > 2) is considered. 

Theorem 2.2. A sujfieient eondition for © is 

det{{F'^)-^nK{F'^)-T)\ 
or 


< 


( 6 ) 


det 

Proof For i ^ j, the left-hand side of @ can be rewritten as 


< min 1 + —, 
d 


^-2 



(7) 


qJiF'KFBMr^q 


^ det{{F^)-^BK{F^)-^)-^I + {F'k)-^Bk{F'k)-^ - det((q,)-iDi^(q,)-^)3/ 


1 


<1 ^-T^4 


= - ^ det((Fi,)-iD^(Fi,)-'^)3 + qj [{F'j,)-^Bk{F'j,)-^ - deti{F^)-^BK{F'j,)-^pI 


Qi, 
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where we have used Qj = —\ for the equilateral simplex K. Prom this, we have 


<_ldet((F;,)-iB,,(F;,)-^)U||(F^)-'B^(P^)-"^-det((P^)-iB^(P^)-^)^/||, 


where || • || is the matrix 2-norm. Thus, a sufficient condition for 0 i 


IS 


-- det((F^)-iDK(F^)-"’)3 + \\{F'k)-^^k{F'k)-^ - det{{F'K)-^^K{F'K)-^)-^I\\ < 0. 

From this we can obtain 

Denote the eigenvalues of (T^) ^Dk(F^) ^ by Ai > • • • > > 0. Then condition 0 is equivalent 

to 


Ai 


(Al • • • \d)d 

Meanwhile, 0 is equivalent to 


< min 




1 

'Al 



■f 


(Al • • • Xd)d 


^d_ 


1 



1 

< - 
- d 


Xi 


(Al • • • Ad)d 


r-1 


< - i = ,d 
d 


1 Xd 

1-u <- - —r < - 

^ iXi---Xd)d iXi---Xd)d 


Al , 1 


Thus, 0 , or 0, implies the right inequality of Q. 

To prove the left inequality of 0 , we start from ([^ and have 


l\-d^ 


Al 




Al _ / Al \ d 

17^“ Vw ’ 


(Ai---Ad)d xdx^d ^Arf. 


or 


> ( 

1 - b 

\XiJ - V 

dJ 


Using this, we have 


Xd ^ Xd 


^ _ {'^d\ 

h - ~ \xJ 


(Ai---Ad)S Xjx^ 
Thus, 0 implies the left inequality of Q. 

Notice that the bound in Q has the value 


A<1\^ >1_ 1 
“ d 


for 2D {d = 2) 
.225, for 3D (d = 3). 


( 8 ) 


(9) 


□ 
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Moreover, ^ is sufficient for Q because ^ implies that all of the eigenvalues of ^ 

are equal to ck and thus the left-hand side of 0 is equal to one which is less than the right-hand 
side. This means that Q is weaker than ([5). 

Unfortunately, Q is still stronger than (^. Consider the triangular and tetrahedral elements in 
Fig. [^and the case with D = /. It is easy to see that the elements are non-obtuse and satisfy A 
direct calculation shows that 


1.732, for 2D elements in Fig. [^a) 

2.151, for 3D elements in Fig. [^b) 

which violates 0. Nevertheless, as will be seen in the next section, the conditions 0 and 0 provide 
useful guidelines for generating meshes that improve the MP satisfaction of the scheme. 


max ■ 




\-T\ 


KJ 





3 Anisotropic mesh adaptation 


In this section we study anisotropic mesh adaptation for the anisotropic diffusion problem 0 and 
0. We use the so-called M-uniform mesh approach 123130] where an adaptive mesh is viewed as a 
uniform one in the metric specified by a tensor M = M(x) which is assumed to be symmetric and 
uniformly positive definite on D. For the moment, we assume that M has been given. Its choice will 
be discussed later. It is shown in 1301 that an M-uniform tetrahedral mesh Th satisfies 


\K\ det(Mi^)^ = ^, ^Ke Th (10) 

itr {iF'K)-\MK)-\F^)-^) = det {{F'j,)-\Mk)-\F'k)-^)K VK G % ( 11 ) 


where 


= ~i^\ f M(a:)da;, 
|A I Jk 


CTfe = det(Mi^)2 

KeTh 


( 12 ) 


Condition (10) is called as the equidistrihution condition which requires all of the elements to have 
the same size in the metric and therefore determines the size of K from det(Mx)2. On the 
other hand, 0 is called the alignment condition which requires measured in the metric M^, 
to be similar to the reference element K that is measured in the Euclidean metric and thus controls 
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the shape and orientation of It can also be shown that the principal axes of the circumscribed 
ellipsoid of K are required to be parallel to the eigenvectors of while their lengths are reciprocally 
proportional to the square roots of the respective eigenvalues m- 

M-uniform or nearly M-uniform meshes can be generated using various strategies; see m Section 4] 
for more detailed discussion. We mention just a few here, the Delaunay triangulation method O [6l O 
[53] . the advancing front method miiiH], the bubble mesh method [Mj, the combination of combining 
refinement, local modification, and local node movement 13 El im EH 05], and the variational method 
1221 . Particularly, we mention two C++ codes which can directly take the user supplied metric tensor. 
One is BAMG (Bidimensional Anisotropic Mesh Generator) for 2D meshes developed by Hecht [25] 
based on the Delaunay triangulation and local node movement. The other is MMG3D (Anisotropic 
Tetrahedral Remesher/Moving Mesh Generation) developed by Dobrzynski and Frey [16] based on 
refinement and local node movement. The latter is used in our computation. 

A key component of the M-uniform mesh approach is to define the metric tensor. We consider four 
choices here. The first one is 

= /, Vx G D. (13) 

This is the simplest choice, and the resulting meshes are uniform or nearly uniform. 

The second choice is based on linear interpolation error. It is known that the error for the linear 
finite element solution to 0 and Q is bounded by the error in the piecewise linear interpolation of 
the exact solution. Thus, it is reasonable to define the metric tensor based on linear interpolation 
error. A metric tensor based on minimization of the semi-norm of linear interpolation error is 
given [29] by 


Madap{K)= I+—\HKiu'^)\ det (I +—\Hk{u^)\ 


1 ^ 
■ 5 


I+—\Hk{u'^)\ 


\JK G Th (14) 


where is the finite element solution, Hk{u^) is a recovered Hessian of over iF, \Hk{u^)\ is the 
eigen-decomposition of Hk{u^) with the eigenvalues being replaced by their absolute values, and 
is defined implicitly through 

^ \K\^det{Madap{K))^2\^\. (15) 

KeTh 

This equation determines ah uniquely and can be solved, for instance, by the bisection method. 
Moreover, the choice of ah in this way concentrates roughly 50% of the mesh points in the region 
where the solution changes significantly. 

The third and fourth choices are related to the diffusion matrix. Recall that 0 requires that all 
of the eigenvalues of the matrix (F^)“^(Mx)“^(F^)“^ be equal to each other. It is mathematically 
equivalent to 

{F'k)-\Mk)-Hf'k)-^ = ckI (16) 

for some constant ck on K. Comparing this with 0 suggests that we choose the metric tensor as 

M{K) = 0Kn-\ \/KeTh (17) 


where 9k is an arbitrary positive piecewise constant function and is the average of D over K. From 
(16), it is not difficult to see that any M-uniform mesh associated with this metric tensor satisfies 0 
and therefore, the finite element solution satisfies MP. 

Then, the third choice is 


Mdmp{K) ^KeTh 


(18) 
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which corresponds to = 1. For the fourth choice, we take advantage of the arbitrariness of 9k in 
0 and choose it to minimize a bound of semi-norm of linear interpolation error. This way we 
can combine the desire of preserving MP with mesh adaptation. The metric tensor reads as (cf. [40] ) 

2 

^DMP+adap{K) = (l-\ - Bk] det (Dx) 3 (19) 

V J 


where 


Bk = det (Bj^)-5 ||B-i|| • \\B)K\HK{u^)\f, 


The M-uniform meshes associated with M.dmp and M^oMPpadap satisfy the alignment condition 0, 
and the mesh elements are aligned along the principal diffusion direction defined as the direction of 
the eigenvector corresponding to the largest eigenvalue of the diffusion matrix D. 

It should be pointed out that in practical computation, it is difficult, if not impossible, to generate 
a perfect M-uniform mesh that satisfy the conditions (10) and 0 . Thus, it makes sense to measure 
how far a given mesh is from satisfying these conditions. From 0 and ( [IT| ), we can define the 
equidistribution and alignment measures as 


Qe^{K) 


Qali{K) 


A^|ii:|det(Mi^)^ 


( 20 ) 

( 21 ) 


It is not difficult to show that the maximum norm of both functions, ||Qeg||L<^ and \\Qali\\L°^^ has 

the range [l,oo). Moreover, \\Qeq\\ poo 1 and WQaiiW poo 1 imply a perfect M-uniform mesh. The 

larger ||Qeg||L'^ and ^he farther the mesh is away from being M-uniform. It is worth 

emphasizing that only Qau is related to MP satisfaction. The quantity Qeq indicates how evenly the 
error is distributed among the elements. 

It is interesting to point out that the definition ( [2l| ) is slightly different from a more common 
definition m where the trace of (F^) ^{Mk) is used. 


4.(k)^ ^tr((^^)-nM.)-nF;.)-q _ 

det((Fj,)-i(M,,)-i(F^)-q^ 


( 22 ) 


These two definitions are equivalent since the trace and 2-norm of a positive definite matrix are 
equivalent. The main motivation for which we use (21) is that the condition ([^ can now be rewritten 
as 

1\-^' 


IICaHlI 


L°° 




(23) 


As mentioned before, this condition is hard to satisfy. Indeed, for meshes shown in Fig. [T} we have 
IlCaiill poo - 1.732 (2D) and 2.151 (3D). Although they violate ([^, the meshes can still be considered 


to be close to satisfying the alignment condition since is relatively small. Moreover, (23) 

shows a connection between the alignment condition and the preservation of MP. 
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4 Numerical examples 


In this section we present two three-dimensional examples to demonstrate the performance of the 
anisotropic mesh adaptation strategy described in the previous section with the four metric tensors. 
An iterative procedure for solving PDEs using anisotropic adaptive meshes is shown in Fig. The 
PDF is first solved on the current mesh using piecewise linear finite elements. Then, the metric 
tensor M is computed, followed by the generation of a new mesh for the metric tensor using MMG3D. 
This procedure is repeated a few times untill the quliaty measures (20) and (21) cannot not be 
improved further. In our computations, the procedure is repeated five times and numerical experiment 
shows that there is no significant improvement in the results when more iterations are used. The 
computations are performed in the framework of an open source finite element software, FreeFem++ 
developed by Hecht [26], with the linear solver being chosen as the conjugate gradient method with 
tolerance 10“^^. Moreover, the Hessian of the finite element solution used in the computation of 
the metric tensors is computed using the “mshmet” library embedded in FreeFem++. The “mshmet” 
library first approximates the gradient (first derivatives) using the least squares fitting, then uses the 
approximated gradient in the least squares fitting for the Hessian matrix (second derivatives). 






Recover solution 


Generate 

Given a mesh 

-> 

Solve PDF 

-> 

derivatives and 

-> 

new mesh 





compute M 


according to M 


Figure 2: An iterative procedure for adaptive mesh solution of PDEs. 


The meshes associated with the four metric tensors will be denoted with the same notation for the 
metric tensors. For instance, a mesh associated with will be called an mesh. We consider 

the diffusion matrix in the form 



sin (j) cos 9 

— sin 9 

cos (j) cos 9 


'ki 

0 

o' 


sin (j) cos 9 

— sin 9 

cos (j) cos 9 

B = 

sin (j) sin 9 

cos 9 

cos (j) sin 9 


0 

k2 

0 


sin (j) sin 9 

cos 9 

cos (j) sin 9 


coscj) 

0 

— sine/) 


0 

0 

k3 


coscj) 

0 

— sine/) 


where ki is the dominant eigenvalue, (j) is the angle between the principal diffusion direction and 
the positive z-axis, and 9 is the angle between the projection of the principal diffusion vector on the 
x^-plane and the positive x-axis. 


Example 4.1. For the first example, we would like to compare the L^-norm of solution errors 
obtained from different meshes and check the violation/satisfaction of MP. We consider problem ^ 
in the unit cube D = (0,1)^. ID is chosen in the form of (24) with cj) = —7r/4, 9 — Stt/G, ki — 100, 
k 2 = 10, and /cs = 1. The source function / and the boundary function g are chosen such that the 
exact solution is given by 

u = e-ioo((^-0-5)"+(y-o.5)2-o.i2) ^ 2 ^ ^25) 


The continuous problem satisfies MP and the solution is in the interval (0,1 + e]. 

The exact solution in the domain D and on the cross-sections x — 0.5 and 2 : = 0.5 is shown in 
Fig. [3} Fig. 13 shows the four types of meshes with a view of the inner mesh by cutting a corner of 
the domain. Fig. shows the meshes on the cross-section x — 0.5. The elements in the M^bmp niesh 
are aligned well along the primary diffusion direction 9 — 57r/6, while the elements in the mesh 





















are aligned along the direction of 6 = 7r/4. The elements in the M^adap niesh concentrate around the 
central region where the solution changes rapidly. The elements in the M.jjMP^adap mesh not only are 
aligned along the primary diffusion direction but also concentrate around the central region, which is 
a result of combining MP preservation and adaptivity. 

Table shows the L^-norm of the solution error obtained on Madap^ ^dmp^ and M^MPpadap 
meshes. The minimal value in the computed solution, which is denoted by Umin and indicates the 
undershoots and violation of MP, is also reported in the table. The convergence history of the L^-norm 
of the error is shown in Fig. It is clear that the error converges at a rate of second order as the 
mesh is refined. Moreover, the numerical solution obtained using the mesh has the smallest 

error, which is consistent with the fact that Madap is formulated based on the minimization of the 
interpolation error (and the finite element error). The solution obtained using the 'Mjjmp mesh has 
the largest error due to the fact that the mesh is designed to satisfy MP, which is generally in conflict 
with error reduction. The error of the solution obtained using M^jjMPpadap is between those using 
Madap 2 tnd M^dmP: which is expected from the construction of M^jjMP+adap ^md shows a good balance 
between mesh adaptivity and MP preservation. 

It can also be observed from Table that the numerical solution obtained from and Madap 
violate MP even for very fine meshes. On the other hand, the numerical solutions obtained from 
^DMP and M^MP+adap violates MP for coarse meshes but satisfy MP when the mesh is sufficiently 
fine. This indicates that 'M^mp and M^^MP+adap can help improving the satisfaction of MP for the 
finite element solution. 



Figure 3: Example 


4.1 


Exact solution in the domain ft and on the cross-sections y 


0.5 and z = 0.5. 


Example 4.2. For this example, we would like to discuss the relation between mesh quality 
and MP satisfaction for three different cases. The domain is (0,1)^ with a cubic hole [0.4, 0.6]^ cut 
inside. The problem is in the form 0 with / = 0 and the Dirichlet boundary condition = 0 on 
the outer surface and u — A on the inner surface. The diffusion matrix is in the form of (24) with 
(p — 7T/ki = 100, k 2 — 10, and k^ — 10. The primary diffusion direction is on the plane oi y — z = d 
(—1 < d < 1) and in the direction defined by the angle 9. 

We first consider two different values of 9 and then a diffusion matrix with a strong anisotropic 
feature. The continuous problem satisfies MP and the exact solution satisfies 0 < < 4. 


Case 1. d = 7t/A. Fig.j^shows the numerical solution obtained using a fine mesh of 5,952,000 
tetrahedra. Four different meshes are shown in Fig. Table shows the minimal value Umin in 
the numerical solution obtained using different meshes. The corresponding mesh quality measures 


9 









(a): mesh, N — 24,576 


(b): Madap mesh, N = 17,304 



(c): M£)mp mesh, N — 29,023 



(d): MoMP+adap mesh, N — 20,099 


Figure 4: Example 


4.1 


Different meshes with a view of the inner mesh by cutting a corner. 


IIQa/illiZ, WQaliW L/OO ^ and ||Qe(7llL2 are also reported in Table Recall that only Qau is related to 
MP satisfaction, while Qeq indicates how closely the mesh satisfies the equidistribution condition (10). 
The mesh quality measures for mesh are calculated using M = in order to see if it satisfies 
(23). The numerical solutions from all the meshes have the maximal value of Umax = 4. 

It can be seen that the numerical solution obtained using the Madap mesh has undershoots {umin <0) 
while the solutions obtained using other meshes satisfy MP. Interestingly, for this case MMG3D 
produces the same mesh from 'Momp and An explanation is that the initial mesh, the mesh, 
is already very close to an M^mp mesh (with ||(5a/7||L2 = 2.32). Since a half of the mesh elements 
are aligned with the primary diffusion direction and the anisotropy is not significant, the elements 
do not need to be stretched more. In fact, the mesh is closer to an Mdmp mesh when ki = 50 
with = 1.72. When the anisotropy is more significant, the mesh is farther away from 

the M.JJMP mesh and the elements need to be stretched more along the primary diffusion direction. 
MMG3D will then adapt the mesh to be more different from the mesh, as will be shown in Case 3. 
Moreover, both Madap ^md M^jjMP+adap meshes have elements of worse shape (with large ||Qa/z||L°°) 
but the overall quality is acceptable (with relatively small ||Qa/z||L2)- 
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(c): 'Momp mesh 


(d): MoMP+adap mesh 


Figure 5: Example 


4.1 


Different meshes at cross-section x = 0.5. 


Case 2. 9 = Sir/A. In this case, the elements of the mesh are no longer aligned along the 
primary diffusion direction. Meshes on the cross-section of z = 7/ are shown in Fig. Notice that 
the orientation of elements in the M^jjmp mesh (Sn/A direction) is very different from that in the 
mesh {tt/A direction). The results of Umin and mesh quality measures are listed in Table 

One can see that none of the meshes satisfies (23), however, recall that ([^ or 0 is only a sufficient 


condition for the MP satisfaction. Notice that the numerical solutions obtained from both the Mjjmp 
mesh and M^MP+adap mesh satisfy MP whereas others do not. Fig. also shows that MoMP+adap 
mesh not only tempts to make elements to be aligned with the primary diffusion direction but also 
concentrates elements to minimize the solution error. 


Case 3. 9 — tt/A^ ki — 1000, ^2 = 1, fcs = 1. In this case, the diffusion is much faster in the 


direction of 0 = 7r/4 than in other directions. Fig. 10 shows the numerical solution using a fine 
mesh with 5,952,000 tetrahedra. A planar view of four different meshes at cross-section z = y is shown 
in Fig. IT The results of Umin and mesh quality measures are reported in Table 

Overall, all but the mesh are close to being M-uniform, with relatively small ||Qa/z||L2 and 
||(5eg||L2- But they do not satisfy ^ so there is no guarantee that the solutions will satisfy MP. 
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Table 1: Numerical results obtained with different meshes for Example 


4.1 



N 

384 

3,072 

24,576 

196,608 

750,000 

1,572,864 

mesh 

L? error 

1.99e-l 

1.07e-l 

5.86e-2 

1.98e-2 

8.85e-3 

5.55e-3 


'^min 

-2.93e-2 

-4.06e-2 

-2.02e-2 

-4.79e-3 

-9.71e-4 

-3.87e-4 

^^adap 

N 

572 

2,865 

17,304 

133,012 

627,329 

1,521,648 

mesh 

P error 

1.20e-l 

5.67e-2 

1.41e-2 

3.72e-3 

1.53e-3 

9.56e-4 


'^min 

-1.94e-2 

-2.26e-2 

0 

-4.17e-4 

-5.68e-5 

-7.04e-5 

Mdmp 

N 

1,254 

5,962 

29,023 

201,039 

736,646 

1,542,910 

mesh 

P error 

2.74e-l 

1.96e-l 

9.84e-2 

3.40e-2 

1.59e-2 

1.03e-2 


'^min 

-2.90e-l 

-1.28e-l 

-7.84e-2 

-2.52e-3 

0 

0 

MoMP+adap 

N 

1,017 

3,500 

20,099 

132,215 

477,802 

978,845 

mesh 

P error 

2.71e-l 

1.48e-l 

5.12e-2 

1.45e-2 

6.34e-3 

3.88e-3 


'^min 

-6.23e-2 

-2.25e-2 

-5.60e-2 

0 

0 

0 


Indeed, the numerical solutions obtained from all four meshes violate MP except for fine Mjjmp and 
^DMPpadap mcshcs. As Can be seen from Table M^omp and MjjMP+adap meshes improve the 
MP satisfaction of the numerical solution. For the numerical solution obtained using the mesh 
with N = 380,928, it has the minimal value Umin — —7.63 x 10“^, and for the 'Madap mesh with 
N = 315,947, Ujnin — —5.79 x 10“^. On the other hand, the numerical solutions obtained using the 
Mdmp mesh with N = 310,147 and the mesh with N = 396, 625 already satisfy MP. 


5 Application to fractured reservoir simulation in petroleum engineering 

Numerical simulation plays an important role in petroleum engineering to predict production rate, 
optimize hydraulic fracturing design, and evaluate enhanced oil recovery processes. In those computa¬ 
tions, the mesh has to be sufficiently refined around wellbore or fractures to accurately represent the 
flow effects thereon because the fractures are significantly smaller but with much higher permeability 
comparing to the reservoir matrix. It is challenging to consider the mesh inside a fracture due to its 
small scale. 

In recent years, shale gas reservoirs have become an attractive source for natural gas production 
largely due to the advancement of horizontal drilling and hydraulic fracturing techniques [60] . Shale 
reservoirs typically have extremely low permeability. For example, the representative permeability is 
1.0 X 10“^ millidarcies (mD) in the Barnett shale but is 5.0 x 10^ mD in fractures [52]. The high 
permeability in fractures makes it possible to produce the gas from the shale while it makes the 
reservoir highly heterogeneous and anisotropic. The complexity of fracture network together with the 
complexity of shale pore structure makes shale-gas reservoir simulation a very challenging task. Some 
researchers focus on fluid flow models in nano-scale shale pores [m Ha sa EH and some other focus 
on different models of fracture network in the reservoir [IIlll9l|52l[56]. 

As our first exploration of the use of anisotropic mesh adaptation in fractured reservoir simulation, 
we consider the steady-state fluid flow that can be applied to effective permeability calculation [41159]. 
For simplicity, we consider incompressible single-phase fluid flow and choose the single porosity dual 
permeability model. Darcy’s Law is chosen to describe the flow in both matrix and fractures with 
different effective permeability in the corresponding regions. We also assume that permeability does 
not change with pressure. Under this setting, we can perform mesh adaptation in both the matrix 
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N: number of elements 


Figure 6: Example 4.1 Convergence history of L^-norm of the error for different meshes. 


and the fractures. 

We consider a reservoir with 3000 ft in x-direction, 1000 ft in ^-direction, and 300 ft in 2 ;-direction, 
and denote the domain as fl. Fig. 12 shows the sketch of the half-panel of the reservoir {0 < y < 500). 
A horizontal well is located in the center of the reservoir with length — 1400 ft along x-direction. 
The radius of the wellbore is 0.3 ft. The reservoir pressure is — 3800 psi and the pressure in 
the wellbore is = 1000 psi. There are three fractures with half-length Lf — 400 ft and width 
Wf = 0.01 ft at the location x = 800 ft, x = 1500 ft, and x = 2200 ft, respectively. The angle 
between the fractures and the positive x-axis is 7 r/ 4 . The height of the fractures is the same as the 
height of the reservoir. The subdomain formed by the fractures is denoted as flj. The permeability 
in the matrix is kmpar = 0.1 mD in the direction parallel to the fractures, is kmperp = 5 mD in the 
direction perpendicular to the the fractures on x^-plane, and is kmpar in fhe direction of z-axis. The 
permeability in the fractures is kpf = 10^ mD along the fractures and is considered the same as that 
permeability in the matrix in other directions. For successful computation, a common technique is to 
scale up the width of the fracture while keeping the same conductivity, that is, KpfWf = 100 mD-ft 
for this example. In our computations, the width of each fracture is scaled up to Wf — 10 ft and the 
permeability is reduced to Kpf = 10 mD. Outside of the fractures, the primary diffusion direction in 
the matrix is perpendicular to the fractures. 

The faces on the left side (x = 0), right side {x = 3000), and back side {y = 500) are denoted as 
Fi. F 2 denotes the front side = 0), and F 3 denotes the bottom side {z = 0) and top side {z = 300). 
F 4 denotes the faces of the fractures on the front side F 2 . The horizontal well is in the center of the 
face F 2 along the x-direction. For simplicity, we assume that the pressure on the whole surface F 4 is 
the same as the pressure in the wellbore Pyj — 1000 psi. The pressure on Fi is always the same as the 
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Figure 7: Example 


4.2 


Case 


1. Numerical solution in the domain ft and on the cross-section z = y. 


reservoir pressure Pr — 3800 psi. There is no flow through r 2 except in r 4 , and no flow through Fs 
either. The mathematical formulation of the model is given by 


with 


-V- (KVF) = 0, 

in D 

P = Pr, 

on Fi 

P^Pw, 

on F 4 

II 

on (F 2 \F 4 ) U F 3 



■7.5 

2.5 

0 ■ 



■ 2.55 

-2.45 

0 ■ 

]K = 

2.5 

7.5 

0 

in Dj, 

and K = 

-2.45 

2.55 

0 


0 

0 

0.1 


0 

0 

0.1 


in 


( 26 ) 


( 27 ) 


where the viscosity of the fluid and the porosity of the matrix are taken out from the equation from 
the assumption that they are independent of pressure. Gravitational effects are also ignored. The 
numerical solution with a fine mesh of 2,804,175 tetrahedra is shown in Fig. 13, where the mesh is 
manually refined around the fractures. 

The xy-planar view of four different meshes ^dmp^ and M.j)MP+adap at cross-section 

z = 100 ft are shown in Fig.[^ In this example, the mesh is the initial mesh that has manual local 
refinement around the fractures. The numerical solutions obtained from and M^adap meshes have 
pressure values larger than 3800 psi inside the domain, which violates MP and causes spurious back 
flow as shown in Fig. ^ The red shaded regions indicate unphysical pressure values. The numerical 
solutions obtained using M^mp and M.p)MP+adap meshes both satisfy MP, and the gradients of pressure 
are shown in Fig. 

Different meshes at cross-sections along the fracture at x = 1500 ft and around the fractures at 
cross-section of 2 ; = 100 ft are shown in Figs. 17 and 18, respectively. 'Mdmp and 'MjjMPpadap 


meshes clearly are aligned along the principal diffusion direction both in fracture and in the matrix 
(where the primary diffusion direction is perpendicular to fractures). The alignment of mesh elements 
helps balance the anisotropic diffusion, and the numerical solutions obtained using both the M^mp 
and M^oMP+adap meshes satisfy MP. M.]jMP+adap mesh also shows a degree of element concentration 
around the fractures. Table lists the maximum pressure values for different meshes, which shows 
that improper meshes can lead to unphysical solutions in the fractured reservoir simulation. 
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(a): mesh, N — 47,616 
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(b): Madap mesh, N — 36,066 



(c): Momp mesh, N = 47,616 (d): MoMP+adap mesh, N — 35,528 


Figure 8 : Example 4.2 Case 1 . Different meshes at cross-section z = y ior 9 — tt/A. 


6 Conclusions and comments 


In the previous sections we have studied anisotropic mesh adaptation and maximum principle preser¬ 
vation for the finite element solution of three-dimensional anisotropic diffusion problems. We have 


defined a quality measure Qaii in (21) that provides a good guidance for MP satisfaction. And we 


have also developed some new sufficient conditions (cf. (©, 0, and d^) that guarantee the linear 
finite element solution to satisfy MP. Comparison has been performed for four different metric tensors 
that are used to generate adaptive meshes with the existing software MMG3D. 

One of the metric tensors is the identity matrix which typically leads to uniform or almost uniform 
meshes. The other three are related to the diffusion matrix and/or the finite element solution. 

(14) is based on minimizing semi-norm of linear interpolation error. M^adap meshes concentrate ele¬ 


ments near the regions where the Hessian of the solution is large, and gives the smallest error. M^mp 


(18) is taken as the inverse of the diffusion matrix and leads to meshes with elements aligned along 
the primary diffusion directions, which improves the MP satisfaction of the finite element solution but 
gives the largest error among all of the considered metric tensors. The last metric tensor, M^j^MP+adap 


(19) is proportional to the inverse of the diffusion matrix, with the coefficient of proportionality being 
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4.2 


Table 2: Minimal solution values and mesh quality measures using different meshes for Example 
Case 1. 


Mesh 

N 

"^min 

1 \Qali 1 II^ 

\\Qali\\L°° 

1 iQegl | l 2 


744 

0 

2.32 

2.58 

1.00 

Mjd 

5,952 

0 

2.32 

2.58 

1.00 

mesh 

47,616 

0 

2.32 

2.58 

1.00 


380,928 

0 

2.32 

2.58 

1.00 


3,047,424 

0 

2.32 

2.58 

1.00 


976 

0 

3.22 

13.0 

1.15 

^^adap 

5,616 

-1.47e-4 

3.98 

22.17 

1.06 

mesh 

36,066 

-2.38e-3 

4.38 

47.11 

1.11 


278,023 

-1.30e-3 

4.25 

106.7 

1.16 


2,779,767 

0 

4.21 

151.5 

1.38 


744 

0 

2.32 

2.58 

1.00 

M.DMP 

5,952 

0 

2.32 

2.58 

1.00 

mesh 

47,616 

0 

2.32 

2.58 

1.00 


380,928 

0 

2.32 

2.58 

1.00 


3,047,424 

0 

2.32 

2.58 

1.00 


748 

0 

2.28 

3.36 

1.08 

^^DMPPadap 

6,539 

0 

2.32 

7.41 

1.07 

mesh 

35,528 

0 

2.27 

9.01 

1.17 


397,125 

0 

2.29 

8.14 

1.12 


3,361,403 

0 

2.32 

12.43 

1.05 


taken as a function to minimize the semi-norm of linear interpolation error. It provides a good 

balance between MP preservation and mesh adaptivity (Madap)- Meshes associated with 

this metric tensor improves the MP satisfaction of the finite element solution while keeping the error 
minimal. 

The anisotropic mesh adaptation procedure is applied to fractured reservoir simulation in petroleum 
engineering. Numerical results show that mesh (with manual adaptation around fractures) can 
lead to unphysical solutions, Madap is capable of concentrating mesh elements around the fractures 
to reduce solution errors but may still violate MP, while both 'Mjjmp and M^^MP+adap tend to align 
elements along the primary diffusion direction where the permeability is significantly larger than those 
in other directions and produce no artificial oscillations in the computed solution. 

Overall, the numerical observations we made here for three dimensional problems are consistent 
with those for two dimensions reported in m- However, numerical experience suggests that it is 
much harder in 3D to make the mesh to satisfy or closely satisfy the alignment condition (0 . This is 
especially true for the mesh adaptation situation (with the metric tensor depending on the solution) 
for which QaU is relatively large for some elements although (which is an average of Qau) 

stays relatively small. How to generate better nearly M-uniform meshes for a given metric tensor in 3D 
has attracted interest from researchers (e.g. see HSIIIH]) and certainly deserves more investigations. 


16 

















(a): mesh, N = 47,616 


(b): Madap mesh, N = 37,899 




N 


X 



(c): Mdmp mesh, N = 47,454 


(d): MDMPpadap mesh, N = 48693 


Figure 9: Example 4.2 Case 2. Planar view of meshes at cross-section z = y for 9 = 37r/4. Notice that 


a planar cut of a 3D mesh does not necessarily form a 2D mesh. 
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4.2 


Table 3: Minimal solution values and mesh quality measures using different meshes for Example 
Case 2. 


Mesh 

N 

"^min 

1 \Qali 1 II^ 

\\Qali\\L°° 

1 iQegl |l2 


744 

0 

7.15 

9.42 

1.00 

Mjd 
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7.15 
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mesh 
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7.15 

9.42 
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-1.16e-5 

7.15 

9.42 

1.00 


3,047,424 

-5.10e-7 

7.15 

9.42 

1.00 


834 

-1.26e-l 

3.87 

22.66 

1.24 

^^adap 

5,485 

-3.10e-3 

4.99 

42.44 

1.34 

mesh 

37,899 

-2.46e-3 

4.42 

35.13 

1.29 


314,228 

-7.62e-4 

12.1 

291.4 

1.10 


2,659,107 

-8.26e-5 

13.2 

491.3 

1.18 


834 

0 

3.46 

14.5 

1.08 

Mdmp 

6,052 

0 

2.43 

15.08 

1.03 

mesh 

47,454 

0 

2.08 

15.3 

1.02 


378,787 

0 

1.93 

22.2 

1.01 


3,035,251 

0 

1.87 

24.0 

1.00 


744 

0 

3.46 

9.42 

1.10 

^^DMP+adap 

6,284 

0 

2.58 

9.42 

1.16 

mesh 

48,693 
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1.10 
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Figure 11: Example 4.2 Case 3. Planar view of meshes at cross-section z = y for 0 = 7r/4, ki = 1000, 


^2 = 1, and /c 3 = 1. Notice that a planar cut of a 3D mesh does not necessarily form a 2D 
mesh. 
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4.2 


Table 4: Minimal solution values and mesh qualities using different meshes for Example 


Mesh 

N 


1 \Qali 1 |l2 

WQaliWh^ 

1 iQef^l |l2 


744 

-6.15e-2 

46.04 

49.99 

1.00 

Mid 

5,952 

-1.09e-l 

46.04 

49.99 

1.00 

mesh 

47,616 

-1.14e-l 

46.04 

49.99 

1.00 


380,928 

-7.63e-2 

46.04 

49.99 

1.00 


3,047,424 

-3.78e-2 

46.04 

49.99 

1.00 


924 

-8.04e-2 

3.72 

23.43 

1.29 

^^adap 

5,137 

-3.38e-2 

4.37 

24.22 
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46,334 
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5.00 

50.52 

1.32 
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73.36 
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94.55 
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3.86 
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1.38 
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0 

2.56 
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1.19 
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0 

2.42 
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1.20 
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0 

10.34 

49.99 

1.74 

^DMP+adap 

6,919 

-2.67e-3 

6.90 

54.49 
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48,267 

-l.Ole-7 

3.77 

63.01 

1.28 


396,625 

0 

2.78 

78.36 

1.22 


2,716,462 

0 

2.43 

54.49 

1.48 


Case 3. 


Table 5: Fractured reservoir simulation. Maximum pressure values obtained with different meshes. 
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N 

P 

max 
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Figure 12: Fractured reservoir simulation. The sketch of the half-panel of the reservoir where the 
angle between the fractures and the x-axis is 7r/4. Fi consists of left side (x = 0), right 
side {x = 3000), and back side {y — 500). F2 consists of front side = 0). F3 consists of 
the bottom side {z = 0) and top side {z = 300). F4 consists of the faces of the fractures on 
front side. 
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Figure 13: Fractured reservoir simulation. Numerical solution in the domain Q and on the cross- 
sections along the fractures at x — 800, x = 1500, and x = 2200. 
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(a): Mid mesh, Pmax = 3821.57 



(b): Closer view of (a) at x = 0 ft 




(d): Closer view of (c) at x = 0 ft 


(c): Madap mesh, Pmax = 3800.07 

Figure 15: Fractured reservoir simulation. Pressure gradients obtained using and Madap meshes. 

The red shaded region have unphysical pressure values that are larger than the reservoir 
pressure. 
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(a): 

Figure 16: 


P-\-adap ruesll, Pmax — 3800 (t^)* P-\-adap ruesll, Pmax — 3800 

Fractured reservoir simulation. Pressure gradients obtained using WI^mp and M^^MPpadap 
mesh. No unphysical pressure values are observed in the computed solution. 
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DMP 


mesh 


(d): MoMP+adap mesh 


Figure 17: Fractured reservoir simulation. Different meshes at cross-section along the fracture at 
a: = 1500 ft. 
























































































































































































































































































































































































